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Abstract 

We present detailed Monte Carlo results for the two-dimensional melting tran- 
sition of various systems up to N = 65536 hard disks. The simulations are per- 
formed in the NVT ensemble, using a new updating scheme. In the isotropic phase 
the bond orientational correlation length £q and the susceptibility %6 ar e measured 
and compared with the predictions of the Kosterlitz-Thouless-Halperin-Nelson- Young 
(KTHNY) theory. From the scaling relation of an d Xe we calculate the critical 
exponent In the phase transition region we use finite-size scaling methods to lo- 
cate the disclination binding transition point and compare the results with the values 
obtained from the behaviour in the isotropic phase. Additionally, we measure the 
topological defect density, the pressure and the distribution of the second moment 
of the local bond orientational order parameter. All results are in good agreement 
with the KTHNY theory, while a first-order phase transition with small correlation 
length and a one-stage continuous transition can be ruled out. 

PACS: 64.70.Dv, 64.60.Fr 

1 Introduction 

The nature of the two-dimensional melting transition is a long unsolved problem 
[|. Melting in two dimensions differs from the three-dimensional case because the two- 
dimensional solid posses only quasi-long-range positional order, while the three-dimensional 
solid is truly long-range positional ordered. This means that the correlation function in 
two dimensions decays algebraically to zero for large distances, while it decays to a non- 
zero value in three dimensions. This absence of conventional long-range order at non-zero 
temperatures in two dimensions was pointed out by Mermin and Wagner ||. Therefore, 
the mean-square displacement of the particles from their ideal lattice position will diverge 
logarithmically with the size of the system and no Bragg peaks in a strict manner can occur 
in the thermodynamic limit. Nevertheless, the other order parameter, which describes the 
bond orientational order, is truly long-ranged ordered, i.e. the orientation of the bonds 
between neighboured particles is correlated over arbitrary distances. 
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There are several theoretical approaches for the description of melting in two dimen- 
sions. Halperin and Nelson as well as Young have developed a theory based on the the 
ideas of Kosterlitz and Thouless Q . The KTHNY theory deals with unbinding scenarios of 
topological defects, where the two order parameters are related to two different topological 
defects: the disclinations and the dislocations. The dislocation unbinding at a temperature 
T m is responsible for the melting transition, while the disclination unbinding at 7] destroys 
the bond orientation. The first continuous transition transforms the solid into a novel hex- 
atic phase, which is short-ranged positional and quasi-long-ranged orientational ordered. 
The second continuous transition transforms this hexatic phase in an isotropic one, i.e. a 
phase with short-ranged positional and orientational order. An alternative scenario has 
been proposed by Chui ||. He presented a theory via spontaneous generation of grain 
boundaries, i.e. collective excitations of dislocations. He found that grain boundaries may 
be generated before the dislocations unbind if the core energy of dislocations is sufficiently 
small, and predicted a first-order transition. This is characterized by a coexistence region 
of the solid and isotropic phase, while no hexatic phase exists. Another proposal was given 
by Glaser and Clark ||. They considered a detailed theory where the transition is handled 
as a condensation of localized, thermally generated geometrical defects and found also a 
first-order transition. Calculations based on the density-functional approach were done by 
Ryzhov and Tareyeva f6j. They derived that the hexatic phase cannot exist in the hard 
disk system. 

Numerical investigations of two-dimensional melting can be done in several ways. On 
the one hand one can simulate the particle system or the defect system 0, on the other 
hand one can study lattice models which describe defects and their elastic interaction || 
or grain boundaries ||. 

The hard disk system is one of the simplest particle models to study the melting tran- 
sition in two dimensions with computer simulation techniques. Even for this simple case 
no consensus about the existence of an hexatic phase has been established. The melt- 
ing transition of the hard disk system was first seen in a computer simulation by Alder 
and Wainwright fllO| |. They used a system of iV = 870 disks, constant volume V and 



molecular dynamics methods (NVE ensemble) and found that this system undergoes a 
first-order phase transition from the solid to the isotropic phase. However, the results of 
such small systems are affected by large finite-size effects. Simulations performed in the 
last years used Monte Carlo (MC) techniques either with constant volume (NVT ensem- 
ble) [JTTJ, [12], [L3|, [14| or constant pressure (NpT ensemble) [|5|, pL6| . Zollweg, Chester and 
Leung [IT]] made detailed investigations of large systems up to 16384 particles, but draw no 
conclusives about the order of the phase transition. The analysis of Zollweg and Chester 



121 for the pressure gave an upper limit for a first-order phase transition, but is compatible 



with all other scenarios. Lee and Strandburg |D| used isobaric MC simulations and found 
a double-peaked structure in the volume distribution. Lee-Kosterlitz scaling led them to 
conclude that the phase transition is of first order. However, the data are not in the scaling 
region, since their largest system contained only 400 particles. MC investigations of the 
bond orientational order parameter via finite-size scaling with the block analysis technique 
of 16384 particle systems were done by Weber, Marx and Binder JT3|. They also favoured 
a first-order phase transition. In contrast to this, Fernandez, Alonso and Stankiewicz JTB[[] 

1 For a critical discussion of this work see Ref. jLTj ]. 
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predicted a one-stage continuous melting transition, i.e. a scenario with a single continuous 
transition at p x = p m and consequently without an hexatic phase. Their conclusions were 
based on the examination of the bond orientational order parameter in very long runs 
of different systems up to 15876 particles and hard- crystalline wall boundary conditions. 



Finally, Mitus, Weber and Marx [14| studied the local structure of a system with 4096 



hard disks. From the linear behaviour of a local order parameter they derived bounds for 
a possible coexistence region. 

In a recent letter we published the first results of simulations of the hard disk 
model in the NVT ensemble with up to 65536 particles to answer the question of the kind 
of the phase transition. We showed that the behaviour of the susceptibility xe and the 
bond orientational correlation length £ 6 in the isotropic phase as well as the value of the 
critical exponent t]q coincide with the predictions of the KTHNY theory. Additionally, 
we performed finite-size scaling (FSS) investigations in the transition region and showed 
that these results are also in agreement with the KTHNY scenario. Here we discuss the 
methods in detail and present additional results for the pressure, the topological defect 
density and the distribution of the second moment of the local bond orientational order 
parameter. All results are compared with the predictions of the KTHNY theory. 



2 Algorithm and measurement 

As mentioned before, we used MC techniques and the NVT ensemble for the simulations of 
the hard disk system. The updating was performed with an improved (non-local) Metropo- 
lis algorithm fll9"| . We consider systems of N = 32 2 , 64 2 , 128 2 and 256 2 hard disks in a 
two-dimensional square box. We find that finite-size effects with these boundary conditions 
are not substantially larger than in a rectangular box with ratio : 2 since no simulations 
in the solid phase were made. This point will be discussed later. The simulations were 
performed on a Silicon Graphics workstation and a CRAY T3E. The CPU time for the 
CRAY was of the order of some month per node, where we have used 7 or 8 nodes. Further 
details are described in Ref. 



Careful attention has been paid to the equilibration of all systems. We controlled that 
the expectation values had stabilized over long time. Additionally, we measured some 
autocorrelation times for smaller systems and estimated the values of larger systems 



(for large correlations lengths) by assuming z ~ 2. We spent at least 10% of the time to 
warm up the system. The measurement frequency was between one measurement per 80 
MC 'sweeps'Q for p = 0.82 and N = 64 2 and one measurement per 5000 MC sweeps for 
p = 0.89 and iV = 256 2 since the measurement is expensive compared to the updating 
steps due to the calculation of neighbours, p is the reduced density since we have set the 
disk diameter equal to one in the whole paper. The number of measurement sweeps for all 
performed simulations is collected in Tab. [1|. The measured observables will be discussed 
in the following. 

2 A sweep for the chain Metropolis algorithm is defined as N trials to move chains of particles p9[. 
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Table V. Number of measurement sweeps that were performed with the chain 
Metropolis algorithm. The acceptance rate was between 50% and 70%. '8 x 
1500' denotes 8 independent data sets with 1.5 x 10 6 sweeps. 



p 






sweeps/10 3 
iV = 64 2 N = 128 2 






N-- 


= 32 2 


N = 256 2 


0.820 






7 x 1 1 






0.830 






7 x 1 900 






0.840 






7 x 1040 






0.850 






7 x 1 i nn 






0.855 






7 x 840 






0.860 


8 x 


1500 


{ x you 






0.865 


8 x 


1500 


7 x 640 


8 x 430 




0.870 


8 x 


1600 


7 x 620 


8 x 470 




0.875 


8 x 


1700 


8 x 500 


8 x 430 




0.880 


8 x 


1500 


8 x 600 


8 x 450 


1900 


0.885 


8 x 


1800 


8 x 1200 


8 x 520 


1900 


0.890 


8 x 


1800 


8 x 1100 


8 x 550 


6 x 1900 


0.895 


8 x 


1800 


8 x 800 


8 x 630 




0.897 


8 x 


1800 


8 x 1500 


8 x 650 




0.898 


8 x 


1500 


8 x 480 


8 x 480 


6 x 710 


0.900 


8 x 


1500 


8 x 640 


8 x 370 




0.905 


8 x 


1500 


8 x 590 


8 x 410 




0.910 


15000 


16500 


8 x 2200 





Bond orientational order parameter and susceptibility 



The orientational order of the two-dimensional hard disk system can be described by the 
(global) bond orientational order parameter ipg. The local value of i/jq for a particle i 
located at r*; = (x, y) is given by 



(1) 



where the sum on j is over the iVj neighbours of this particle and 0jj is the angle between 
the particles i and j and an arbitrary but fixed reference axis. Neighbours are obtained 
by the Voronoi construction ||| . The (global) bond orientational order parameter is then 
defined as 



^6 



N 



iv i=l 



(2) 
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We measured the second and fourth moment of ipg, where the former is related to the 
susceptibility byQ 

Xe = N(^ 2 ). (3) 

Bond orientational correlation length 

The bond orientational correlation function is defined as 

* (f _ o = M0*£3>, (4) 

(p(r)p(r')) 

where 

N 

fair) = ^2ti{r-r i )ip(i ! i (5) 

i=l 

denotes the microscopic bond orientational order-parameter density and 

N 

Pif) = ^2S(f-i\) (6) 

i=i 

is the microscopic particle density. In the isotropic phase the bond orientational correlation 
length ^6 was extracted from the 'zero-momentum' bond orientational correlation function 
ge(x). This is defined as 

g 6 {x - x') = - j J dy dy' g s (f - r') , (7) 

where L denotes the box length. 
In practice we use the definition 



1(1 f L rx+Ax/2 \ * / l i-L rAx/2 

9e(x)~(\— dy dxi/j 6 (f)\ — / dy' dx' i/) 6 (r') 

\\JMk JO Jx-Ax/2 ) \l\ k Jo J-Ax/2 



(8) 



where 

fL rx+Ax/2 

N k = dy dxp{f) , (9) 

JO Jx-Ax/2 
fL pAx/2 

N' k = dy dxp(f) . (10) 

JO J-Ax/2 

Therefore, the distance between two particles in x-direction is not exactly x, but lies 
between x — Ax and x + Ax. Nevertheless, assuming a pure exponential behaviour of the 
correlation function g^x) the integration over Ax causes no error. In the simulations the 
value of Ax was given by the length of a cell of the cell structure, i.e. Ax ~ 2/3, where the 
exact value depends on p and N. g%{x) was fitted with a cosh((L/2 — x)/£ 6 ) in the interval 
a^min < x < L/2, where £ 6 and ^ 6 are related by 

' sinh(-L) . (11) 



2£e V2C 



6, 



3 This definition yields a factor 1 — 2/tt in the thermodynamic limit compared to X6 = N({ip§ 2 ) — (ip6) 2 ) 
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N=64 2 p=0.860 N=64 2 p=0.860 




Figure 1: 'Zero-momentum' bond orientational correlation function ge(x) for 
N = 64 2 particles at p = 0.86 (with arbitrary chosen nomalization). The left 
figure shows the exponential behaviour for large distances, where Ax was about 
0.65. The dotted line is the best fit with a hyperbolic cosine ansatz. The right 
figure illustrates oscillations in ge(x) for small distances (with Ax ~ 0.065). The 
line is a guide to the eye. 



To determine the influence of 'excitations' we compare the results for different minimal 
distances x m i n . The correlations are always dominated by the lowest state of the transfer 
operator 'Hamiltonian', so that it was not necessary to omit points. In Fig. [I] we plot the 
correlation function for the N = 64 2 particle system at p = 0.86 (with arbitrary normaliza- 
tion). The left figure shows the correlation function as obtained from the simulation, i.e. 
with Ax m 2/3. The curve shown is the best fit with a cosh-like behaviour. As one can see 
there are no influences of 'excitations'. Although the fit seems to be consistent with the 
data, there are large deviations. The reason is an oscillating behaviour of g§{x) as shown 
in the right figure, where we have chosen Ax ten times smaller. The same oscillations can 
be seen, if we plot the relative deviations between the data of the first case (Ax ~ 2/3) 
and the hyperbolic cosine fit. This is done in Fig. |2|. Since the oscillating length is about 
1, the curve can be smoothed if one chooses Ax ~ 1. Nevertheless, also with Ax si 2/3 a 
precise determination of the correlation length is possible. Systematic errors coming from 
the oscillations are taken into account. These errors get dominant — compared to our 
statistical errors — for small values of ^6- 

Radial bond orientational correlation function 

In the isotropic phase g§{r) is independent of the angle. Therefore, we use the angle 
averaged quantity 

g 6 {r)^(^ 6 *(Q)^ 6 (r))/g(r) (12) 

for an additional calculation of ^6, where g(r) is the (radial) pair correlation function. The 
radial bond orientational correlation function g 6 (r) was fitted for large distances with an 
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Figure 2: Relative deviations of ge(x) (left picture in Fig. 1) from a hyperbolic 
cosine fit. 



ansatz of the form 



In Fig. [| we plot g&(r) for N = 64 2 hard disks at p = 0.86. The left figure shows the 
oscillating behaviour of ge(r). In order to smooth the curves, ge(r) has been averaged over 
a distance of one. This was done in the right figure, where g§{r) was additionally multiplied 
by r Ve in order to compare the data with an exponential behaviour. 

The values of ^ obtained from ge(r) are affected by larger systematic errors (compared 
to the previous method). The reason is that one has to leave out the points with very 
small and very large distances. The first points have to be omitted since the ansatz is not 
valid in this case, while points with r « L are affected by finite-size effects. In contrast 
to ge(x), where the periodicy just leads to a cosh behaviour, g§{r) has no simple periodic 
behaviour. Therefore, one has to omit the points with large r. Nevertheless, we use the 
radial bond orientational correlation function for a determination of the correlation length. 
In all cases both values of £q coincide within the statistical errors. 

Pressure 

The pressure was calculated from the pair correlation function g(r) in the range 1.0 < r < 
1.2. From 200 bins we extracted the contact value of the pair correlation function by fitting 
the data with a power series of sixth order and extrapolating to g{l). The virial theorem 
relates this value to the pressure by |2l|| 



9e{r) 



r m exp(— r/£ 6 ) . 



(13) 




(14) 
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Figure 3: Radial bond orientational correlation function ge(r). The left figure 
shows the behaviour for small distances. The line is just a guide to the eye. In 
the right figure we plot r V6 ge{r) together with an exponential fit (dotted line). 



where A is the closed-packed area of the system, i.e. A = Ny/3/2. Statistical errors were 
calculated by independent data sets and by performing fits on the whole data sets to a 
Gaussian distribution of g(r) with variance Ag(r). Systematic errors were estimated by 
changing the order of the power series from six to five. 

Our results for the pressure as a function of the system size and the density are collected 
in Tab. |2] and visualized in Fig. |j. The quoted error is the sum of the statistical and 
systematic error. The data show the end of the liquid region and the beginning of a possible 
liquid-solid tie line, while no simulations in the solid phase were made. For densities up 
to p = 0.885, the pressure does not have any finite-size effect within the statistical errors. 
Taking the finite-size dependency of the pressure together with the data of ^6 (which are 
discussed in Sec. Ill), we find that we have reached the thermodynamic limit for the systems 
with N = 128 2 particles up to p = 0.885 and for the systems with N = 256 2 particles up 
to p — 0.89. For densities p > 0.89 there might be still finite-size effects. The results 
are consistent with those of Zollweg and Chester JEfl , who used the same methods but a 
rectangular box with ratio \f?> : 2. Only the value at p = 0.910 shows deviations. This 
could be a result of the square box, which leads to larger finite size effects if the density 
of the system is near the solid phase. Another possibility is that the large systems at 
higher densities are not fully equilibrated. However, this seems to be unlikely due to our 
observation of the pressure as a function of time. Our results give a lower bound for the 
beginning of a coexisting phase of p ~ 0.89, but give neither any conclusive evidence for a 
first-order phase transition or an hexatic phase. It just shows that the compressibility in 
this region is very high. 



Local bond orientational order parameter 



The distribution of the second moment of the local bond orientational order parameter 
| ^6,i | 2 was first studied by Strandburg, Zollweg and Chester [22]. In the case of a first-order 
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Table 2: Pressure for densities in or near the transition region. 



p 




pA /NkT 






N = 32 2 


N = 64 2 


N = 128 2 


N = 256 2 


0.880 


7.803(5) 


7.799(6) 


7.796(7) 


7.795(8) 


0.885 


7.894(5) 


7.900(6) 


7.899(8) 


7.895(9) 


0.890 


7.926(5) 


7.950(7) 


7.950(9) 


7.953(5) 


0.895 


7.910(6) 


7.953(9) 


7.963(9) 




0.897 


7.905(6) 


7.940(6) 


7.956(7) 




0.898 


7.892(6) 


7.934(9) 


7.955(5) 


7.954(4) 


0.900 


7.897(7) 


7.928(9) 


7.951(7) 




0.905 


7.906(6) 


7.901(9) 


7.943(8) 




0.910 


7.916(5) 


7.900(5) 


7.928(5) 






Figure 4: Pressure as a function of the density for various system sizes. Data of 
Ref. [12] are marked by stars. 
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Figure 5: The left figure shows the distribution occurrence of the second mo- 
ment of the local bond orientational order parameter |^6,i| 2 ( m arbitrary units) 
for N = 128 2 hard disks at three different densities. £ indicates the curve, which 
is the linear combination of p\ = 0.89 and p 2 = 0.905. The small inset ampli- 
fies the region with small |^6,i| 2 to show the small difference between the direct 
measurement and the linear combination. Errors are of the order of the distance 
between the two curves. The right figure displays two distributions, which can be 
taken as the initial distributions for the modelling of all others in the transition 
region. 



phase transition (with thin interfaces) one expects that the distribution of the coexistence 
phase is the sum of the fluid, solid and interface distribution weighted with their relative 
areas. On the left picture of Fig. |5] we plot | %/^e,i 1 2 f° r systems with 16384 hard disks at three 
different p's. To check if the distribution at p = 0.898 corresponds to a coexisting phase, 
we compare it with a combination of two other distributions at p\ and p2, respectively. It 
is not necessary that the two chosen densities (pi and p 2 ) are the exact values of lowest 
and highest density of a possible coexisting region. This should work for two arbitrary 
densities, provided that these systems are in the coexisting phase. If a first-order phase 
transitions exists, but the chosen density of p\ = 0.89 is too low or p 2 = 0.905 is too high, 
there are deviations. Obviously, the direct measurement and the modelling are in perfect 
agreement. Moreover, the weights of the two distributions correspond to their theoretical 
values of 8/15 and 7/15, respectively. Nevertheless, an interpretation as the sum of two 
distributions from two different phases of a first-order phase transition makes only sense if 
the system size is larger than the two interfaces. But the results of the following sections 
will show that a first-order phase transition with such small interfaces can be ruled out. 
Therefore the situation is more complicate than in this simple picture. 

The distributions of | ^6,1 1 2 i n the transition region can be modelled as the sum of two 
initial distributions. The reconstruction of these distributions is not unique. A decom- 
position is shown in the right picture of Fig. |. One of the distributions results primary 
from particles with six neighbours, while the other is mainly the sum of distributions from 
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Figure 6: Topological defect density 

^def & 

function of p. Statistical errors 

are too small for a visualization. 



particles with coordination numbers unequal six. 

Additional investigations for the dependency of the distribution of | ip^i | 2 on the system 
size are discussed in Sec. IV. 

Topological defects 

An analysis of the numbers of neighbours of each particle as obtained from the Voronoi 
construction gives a characterization of the defect structure of a two-dimensional system. 
The average number of neighbours is — independent of the state of disorder — six. In a 
perfect solid each particle has six neighbours. Particles with any other number of neigh- 
bours represent a disclination. Dislocations are pairs of disclinations. We define the density 
of defects as 

ndef=^£^\ (15) 

where A^ Nb denotes the number of particles with i neighbours. Alternatively one can take 
the strengths of the disclinations into account and define the density of defects as 

<ef = ^El*-6|Af b = ! E(6 " Nf h . (16) 

i i<6 

However, the difference between the two definitions for p > 0.88 is lower 1%. 

In Fig. P we plot n^ei as a function of p. One can see that there is a linear behaviour of 
ridef for p > 0.890. As in the case of the distribution of iV^il 2 ; this could be explained with 
the coexistence of two different phases. However, if one examines the defect structure of 

4 A linear behaviour of a local order parameter was also found in Ref . Q . 
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several configurations in the transition region and the configurations itself, one finds no hint 
for two coexisting phases, while they are compatible with the picture of a homogeneous 
phase. Therefore, the conventional picture (of a first-order phase transition with thin 
interfaces) of two separated phases is incompatible with the data. The results of Sec. IV 
will confirm this assumption. 

There are different possibilities to explain this linear behaviour. On the one hand there 
could be a weak first-order phase transition with an interface width which is larger than the 
box length L of our largest system of N = 256 2 particles, on the other hand a continuous 
transition with a homogeneous phase. In both cases increasing the density p primary leads 
to a decrease of the defect density, since the average density in the ordered regions in 
higher than the density in unordered regions (i.e. a disclination or dislocation needs more 
space than perfect crystalline structures). In the case of a first-order phase transition the 
defects will form some larger structure, while there is a homogenous distribution for the 
continuous transition. 



3 Simulation in the isotropic phase 

In the isotropic phase we measured the susceptibility and the correlation length of the bond 
orientation. Subsequently, we compare the results with the predictions of the KTHNY 
theory, i.e. a critical exponent 776 of 1/4 and an exponential singularity for the correlation 
length 

Ut)=^exp(b^ 2 ) (17) 

and the susceptibility 

Xe(t) = a x exp(6 x r 1/2 ) (18) 

if t = pi — p — > + . A detailed description of these measurements is given in Ref. [[HJ. 

Our results of x& an d ^6 as a function of the density are summarized in Tab. |3|. We 
analyzed the critical behaviour of Xe{p) an d ^(p) by performing least square fits according 
to Eqs. QT7D and (|Tg). Using all 12 points we got a x 2 P er degree of freedom (d.o.f.) of 0.75 
for £ 6 (t) and 0.65 for Xe(t)^ i- e - the data are in a very good agreement with an exponential 
singularity of the KTHNY type. This is not only a result of large statistical errors as can 
be seen if one uses different approaches for the singularities. For example, a conventional 
second-order behaviour with a power-law singularity of the form ln(£ 6 ) = a — v ln(i) yields 
X 2 /d.o.f.= 4.1. 

All results for the fit parameter are collected in Tab. (|. The values for Xg{p) an d ^(p) 
together with the fitted curves are shown in Fig. |7|. We also made fits where we have 
omitted some data at lower densities. The fit parameters for Xe{p) show only non-essential 
changes, while the changes for ^e(p) are of the order of the statistical errors. An analysis 
of the behaviour of ln(x 6 /£ 6 7//4 ) as a function of ln(£ 6 ) yields the following value of the 
critical exponent: 

77e = 0.251(36) . (19) 
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Table 3: Bond orientational correlation length £ 6 and susceptibility x& f° r various 
densities in the isotropic phase. N refers to the system sizes used. 





p 


N 


£ 6 


X6 


0. 


82 


64 2 


1.513(50) 


3.797(13) 


0. 


83 


64 2 


1.800(35) 


4.693(15) 


0. 


84 


64 2 


2.156(40) 


6.052(24) 


0. 


85 


64 2 


2.635(30) 


8.415(41) 


0. 


855 


64 2 


2.995(35) 


10.30(6) 


0. 


86 


64 2 


3.425(40) 


12.96(9) 


0. 


865 


64 2 , 128 2 


4.14(10) 


17.45(18) 


0. 


87 


128 2 


5.03(15) 


25.00(39) 


0. 


875 


128 2 


6.65(30) 


39.5(8) 


0. 


88 


128 2 , 256 2 


9.56(26) 


75.0(21) 


0. 


885 


128 2 , 256 2 


15.65(51) 


176.8(61) 


0. 


89 


256 2 


38.0(15) 


865(44) 



4 Simulation in the transition region 

We now come to the simulations with p w p ; . Finite-size scaling (FSS) implies for the 
susceptibility 

Xe ~ L 2 - m f(L/£ 6 ) . (20) 

Assuming the prediction of the KTHNY theory, the correlation length diverges at p = p\ 
and / is a constant independent of L. We use this FSS behaviour to locate p ; , where we 
take 7/6 = 1/4. In the hexatic phase (pi < p < p m ) the correlation length £ 6 also diverges, 
so that / is still independent of L. In this phase t/ 6 is a decreasing function of the density 
which goes to zero if p approaches the melting density p m , i.e. at the end of the hexatic 
phase. For p's below p;, one has to take corrections of x& ~ L 2 ~ m for finite correlations 
lengths into account. Our results for the susceptibility are collected in Tab. |5|. 



Table 4: Best fit parameter for the critical behaviour of the correlation length 
and susceptibility. For 0.82 < p < 0.89 we have fitted 12 points, while we used 8 
points in the range 0.855 < p < 0.89. 



Fit 


range 




ln(a) 


b 


Pi 


X 2 /d.o.f. 


Up) 


0.82 < p < 0.1 


39 


-1.44(8) 


0.547(21) 


0.9017(6) 


0.75 


Up) 


0.855 < p < 0J 


S9 


-1.27(13) 


0.505(31) 


0.9006(8) 


0.23 


Xe(p) 


0.82 < p < 0J 


89 


-1.65(3) 


0.847(7) 


0.9002(3) 


0.65 


xa{p) 


0.855 < p < 0, 


39 


-1.60(9) 


0.834(21) 


0.9000(4) 


0.58 
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Figure 7: Susceptibility (full symbols) and bond orientational correlation length 
(open symbols) as a function of density. The curves shown are the best fits 
for a KTHNY behaviour (for all measured points). The critical values of p are 
visualized by vertical lines. 



If we use the FSS behaviour to locate p\ and p m , then the requirement of ^(Pi) = 1/4 
yields 

Pi = 0.899(1) , (21) 

while rj e (p m ) = leads to the estimate p m ^0.91. The value of p x is in agreement with that 
obtained from the singularities of ^{t) and Xe{t)- A slightly different value of T]e (from the 
relation of \& an d ^6 i n the isotropic phase) would not change this situation. Moreover, 
our values of p x are in agreement with the result of Weber et al. |Tt| obtained from the 



fourth-order cumulant intersection (p; = 0.8985(5)). However, it differs from their value 
obtained using the singularity of x& (Pi — 0.913). The result p ; = 0.916(4) of Fernandez et 
al. (TBI is not compatible with our value. 



Another quantity which can be used to analyze the kind of the transition is the fourth- 
order cumulant 

U = 1 - • ( 22 ) 

According to the prediction of the KTHNY theory, U should be independent of the system 
size L in the whole hexatic phase. In contrast to this, in the case of a conventional first- 
order phase transition there is only a single point, where the cumulants of different system 
sizes collapse. Since there is a large region between p x w 0.9 and p m >0.91, the behaviour 
of U can be used to distinguish between a KTHNY and a first-order transition. The 



intersection of the cumulant U in a single point was an argument in Ref. |L3| against the 
existence of an hexatic phase. Unfortunately, statistical errors in our data are too large to 
answer this question as can be seen in Fig. §. 
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Table 5: The susceptibility per particle in the transition region. 



P Xe/N 

N = 32 2 N = 64 2 N = 128 2 N = 256 2 



0.895 
0.897 
0.898 
0.900 
0.905 
0.910 



0.2620(9) 

0.2987(9) 

0.3175(10) 

0.3514(10) 

0.4235(19) 

0.4900(29) 



0.1970(17) 
0.2418(18) 
0.2612(17) 
0.3076(13) 
0.4055(11) 
0.4840(24) 



0.1409(24) 
0.1899(25) 
0.2160(25) 
0.2630(17) 
0.3745(13) 
0.4707(10) 



0.1788(29) 



0.6700 
0.6500 

u 

0.6300 
0.6100 
0.5900 

3.3 3.8 4.3 4.8 5.3 5.8 

L 

Figure 8: Finite-size scaling of the cumulant in the transition region. 
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0.020 



0.015 - 



0.010 - 



0.005 



0.000 




Figure 9: Distribution of the second moment of the local bond orientational 
order parameter (in arbitrary units) for different system sizes. The small inset 
amplifies the region with small | V^s,* 1 2 ? where are the largest deviations. Statistical 
errors are of the order of the symbols in the inset. 



Another possibility to distinguish a first-order phase transition from a continuous tran- 
sition is to study the dependency of the distribution of \4>6 t i\ 2 on the size of the system. 
If the system exhibits a homogeneous hexatic phase, then changing the size of the system 
should not lead to any changes in the distribution. On the other hand, if the transition 
is of first order one would expect that the distribution is a combination of the solid, fluid 
and interface distribution. Therefore, changing the size of the system would result in a 
change of the distribution, because the area of the interface scales only linear with L. In 
Fig. |9| we plot | 2 at p = 0.898 for four different system sizes. Apart from finite-size 
effects, which getting weaker for larger systems, no difference between the distributions can 
be seen. The distributions for the two largest systems coincide within statistical errors. 
Therefore, one can rule out a first-order transition with thin interfaces^, while a first-order 
transition with an interface width larger than the largest system size L and a continuous 
transition are compatible with the data. These results coincide with those of Fernandez 



et al. [16], who performed similar measurements in the NpT ensemble using a rectangular 
box of ratio \/3 : 2. The data of Fig. [| show also that the chosen ratio of the side lengths 
of 1 : 1 causes no large finite-size effects. 



5 The results are also compatible with the occurrence of two very small interfaces (i.e. a width of 0(1)), 
but this can be ruled out due to the examination of the defect structure of several configurations. 



16 



5 Conclusions and outlook 



We presented a detailed Monte Carlo study of the two-dimensional hard disk model in 
the NVT ensemble. The investigations were performed in the isotropic phase and in the 
transition region. 

The behaviour of the defect density as well as the distribution of the local order param- 
eter in the transition region were in good agreement with a simple model of two coexisting 
phases, i.e. the data could be modelled as the sum of two different phases, where the rel- 
ative areas of the two phases are proportional to p. However, the defect structure of the 
system and the distribution of | ip^ 1 2 as a function of L showed, that there are not two 
separated phases with a thin interface. The data can be explained by a weak first-order 
transition with a width of the interface which is larger than the largest system size L or 
by a continuous transition with a homogeneous phase. 

The behaviour of the pressure was compatible with both a first-order and a KTHNY-like 
scenario. The data just give a lower limit of p ~ 0.89 for the coexisting phase. 

In the isotropic phase we examined the dependency of the correlation length and the 
susceptibility on the density p. We showed that the data are in good agreement with the 
prediction of an exponential singularity from the KTHNY theory. The critical exponent 
t]q was derived from the relation on £ 6 and \& f° r ^6 ~~ > °°- We got r/ e = 0.251(36), which 
coincides with the prediction 7/ 6 = 1/4. 

The simulations in the transition region were used to measure the finite-size scaling 
of the susceptibility. The value of p\ = 0.899(1) (assuming r/ 6 = 1/4) coincided with 
those from the KTHNY-like behaviour of ^(p) and Xs{p)- Furthermore, the requirement 
??6(Pm) = led to the estimate p m ^0.91. The data of the fourth-order cumulant U were 
affected by too large statistical errors in order to draw any conclusives. 

In summary, all data are compatible with a KTHNY-like phase transition. A one- 
stage continuous transition (p; = p m ) as proposed in Ref. fl6| and a first-order transition 
with small correlation length can be ruled out[]. Further numerical investigations have to 
be performed to make a clear decision between a weak first-order phase transition and a 
continuous scenario. This could be done for example by studying the positional order in 
the transition region. Work along this line is in progress. 
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